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Recent measurements by the Belle Collaboration of the exclusive production of two charmonia 
in e + e~ annihilation differ substantially from theoretical predictions. We suggest that a significant 
part of the discrepancy can be explained by the process e + e~ — > J/ip + J/ip. Because the J/tp + J/ip 
production process can proceed through fragmentation of two virtual photons into two cc pairs, its 
cross section may be larger than that for J/ip + n c by about a factor 1.8, in spite of a suppression 
factor a 2 /c>s that is associated with the QED and QCD coupling constants. 

PACS numbers: 13.66.Bc, 12.38.Bx, 14.40.Gx 

The charmonium states have long been a key testing ground for our understanding of quantum chromodynamics 
(QCD), both because of their clear experimental signatures and because of the theoretical simplifications that arise 
from their nonrelativistic nature. The huge data sample of continuum e + e _ -annihilation events now being produced 
at the SLAC and KEK B factories makes it possible to investigate charmonium-production processes that have very 
small cross sections. The exclusive production of double-charmonium states is a particularly interesting example 
because the theoretical predictions do not involve any unknown nonperturbative parameters. However, the first such 
measurements by the Belle Collaboration are in substantial disagreement with existing predictions. In this paper, 
we argue that a significant part of the discrepancy between experiment and theory may be attributable to the process 
e + e~ annihilation into J/ip +J/ip, where J/ip is the lowest spin-triplet charmonium state. This process had been 
overlooked because it is suppressed by a 2 /a 2 , where a is the quantum-electrodynamic (QED) coupling, and a s is the 
QCD coupling. However, as we shall show, this process is enhanced by a kinematic factor that is associated with the 
fragmentation of photons into cc pairs. 

The Belle Collaboration has observed e + e _ annihilation into two charmonium states at a center-of-mass energy 
y/s = 10.6 GeV by studying the recoil- momentum spectrum of the J/ip [lj. The collaboration measured the production 
cross section for J/ip + r\ c and also found evidence for J/ip + \c0 and J/ip + r/ c (2S) final states. Recent calculations 
of the production cross section for J/ip + r\ c have given results that are about an order of magnitude smaller than the 
Belle measurement 0, HI • This presents a challenge to our current understanding of charmonium production based 
on NRQCD. 

States consisting of two charmonia with opposite charge conjugation, such as J/ip + r\ c , can be produced at order 
a 2 a 2 through processes e + e~ — > cc + cc that involve only a single virtual photon. (See Refs. 0, Q for examples.) 
States consisting of two charmonia with the same charge conjugation, such as J/ip + J/ip, can be produced at order 
a 4 through the processes shown in the diagrams in Fig. ^ which involve two virtual photons. We find that, in spite 
of its being suppressed by a factor of a 2 /a 2 , the J/ip + J/ip production cross section may be larger than that for 
J/ip + r] c . We suggest that some of the events in Belle's J/tp + rj c signal may actually be J/ip + J/ip events. By taking 
this effect into account, one would decrease the discrepancy between the Belle measurements and the predictions of 
NRQCD. 

The effective field theory NRQCD can be used to write a quarkonium production cross section as a sum of products 
of short-distance coefficients and NRQCD matrix elements Q. The short-distance coefficients can be calculated 
in QCD perturbation theory, but the matrix elements are nonperturbative in nature. In the nonrelativistic limit, 
there is a single independent matrix element for each spin multiplet of charmonium states. The matrix elements 
can be determined phenomenologically from electromagnetic annihilation decay rates. Thus, the cross sections for 
double-charmonium production can be predicted up to relativistic corrections without any unknown nonperturbative 
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FIG. 1: QED diagrams for e + e~ — > 7*7* — > cc + cc. The upper and lower cc pairs evolve into ffi and H2, respectively. 




factors. 

It is convenient to express the cross sections for double-charmonium states in terms of the ratio R, which is defined 
by 



R[H! + H 2 ] 



Ht + H 2 ] 



H+fj.- 



(1) 



where <r[e + e~ — > = 7ra 2 /(3i? 2 cam ), and £?beam = We define the angular variable a; = cos#, where is 

the angle between the J/ip and the beam in the center-of-mass frame. We also introduce a dimensionless kinematic 
variable r = 2m c /£'beam, with m c the charm-quark mass. In the nonrelativistic limit, cross sections for the IS states 
J /if) and r\ c have a common NRQCD matrix element (Oi)is that is defined in Ref. 2\. 

The cross section for J/tp + J/ip receives contributions from the diagrams in Fig. ^ The angular distribution is 



where 



dR 

dx 



[J/1> + J/iA 



16^ 2 a 2 



243m6 
(1 _ r 2)i/2 F ^) 



(Oi)is 



[4(1 - r 2 )(l - x 2 ) + r 



412 ' 



(2) 



F(x) = r 4 x 2 (l~x 2 )[2 + 3r 4 -r b +4x 2 (l-r 4 )} 2 + 2(l + x 2 )(l-x 2 ) 
x [6 - 7r 2 + 4r 4 - r 6 + ix 2 r 2 (l - r 2 )] 2 + 2rV(l - a; 2 ) 
x [3 - 4r 2 + 4r 4 - r 6 + 4xV(l - r 2 )} 2 + 2r 2 [(2 - r 2 ) 2 (3 - 2r 2 + r 4 ) 5 
-ie 2 (72 - 228r 2 + 359r 4 - 292r 6 + 130r 8 - 32r 10 + 3r 12 ) 
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+4a; 4 (9 - 60r 2 + 122r 4 - 112r 6 + 58r 8 - Ur 
+16:r 6 r 2 (l - r 2 )(6 - llr 2 + llr 4 - 2r 6 ) + 64x f V 4 (l - r 2 ) 



„2\2l 



(3) 



The ratio R is obtained by integrating x only from to 1, in order to avoid double-counting of identical final-state 
particles. The angular distribution dR/dx for J/ijj + ip(2S) is given by an expression identical to that in Eq. 
except that one of the factors of (Oi)is is replaced by (Oi/2S, and the range of x is from —1 to 1. 

The cross section for r\ c + r\ c receives contributions only from the two diagrams in Fig.^c) and^d). The angular 
distribution is 



dx 



[Vc + r) c ] 



16n 2 a 2 



(Oi)l s 



243m 6 c 

xr 4 (l-rf¥(l-i 2 ). 



(4) 



The ratio R is obtained by integrating x from to 1. The angular distribution dR/dx for r\ c + i] c (2S) is identical to 
that in Eq. Q, except that one of the factors of (Oi)is is replaced by (Oi)2S, and the range of x is —1 to 1. 



3 



The production of J/ip + r\ c proceeds through e + e~ annihilation into a single virtual photon, which creates a cc 
pair. The second cc pair is then created either by a virtual gluon or by a virtual photon radiated from the first cc 
pair. The cross section was recently calculated by Braaten and Lee 0- The differential ratio is 

+ Vc] = (Oi) lS [3a s r 2 + a(3 + r 2 )] 2 r 2 (l - r 2 ) 3 / 2 (l + x 2 ). (5) 

The ratio R is obtained by integrating over x from —1 to 1. The a 2 term was also calculated recently by Liu, He, and 
Chao 0. It had been calculated previously, but the analytic expression was not given 0. The QED terms increase 
the cross section by about 29% for -Ebcam = 5.3 GeV. 

When -Ebcam is much greater than m c , the relative sizes of the double-charmonium cross sections are governed 
not only by the powers of the coupling constants, but also by the number of kinematic suppression factors r. The 
power-counting rules of perturbative QCD [|J require dR/dx to behave as r 4 in the limit r — > 0, with a further 
suppression factor of r 2 ' Al+A2 if the sum of the charmonium helicities, Ai and A 2 , does not vanish. In the limit r — > 0, 
dR/dx for r\ c + rj c scales as a 2 r 4 , while dR/dx for J/ip + r\ c scales as a 2 r 6 , with the extra factor of r 2 arising from 
helicity suppression. As r — > with x fixed, dR/dx for J/ip + J/ip approaches a constant. The power-counting rules 
of perturbative QCD are evaded because this process has a contribution, corresponding to the diagrams of Figs. Of a) 
and^b), in which two virtual photons, with virtuality of order m c , fragment into two cc pairs. This contribution 
is enhanced because the virtual-photon propagators are of order 1/m 2 instead of order 1/-E< 2 cam . In the amplitude, 
there are also two numerator factors of m c instead of -Ebeam, which arise from the cc electromagnetic currents. Hence, 
the net enhancement of the squared amplitude is (Ehea,m/m c ) . The integrated ratio is further enhanced by a factor 
ln(l/r) because the potential logarithmic divergence in the integral of Eq. J2J as x — » 1 is cut off at 1 — x 2 ~ r 4 . 

Chang, Qiao, and Wang have pointed out that the cross section for e + e~ — * J/ip + 7 is enhanced if -Ebeam is much 
greater than m c 0. The reason for the enhancement is that this process also involves a contribution in which a 
virtual photon fragments into a cc pair. 

We use the results given above to predict the cross sections for producing two S'-wave charmonium states at 
the B factories, with -Ebcam = 5.3 GeV. The ratios R depend on the coupling constants a and a s , the charm- 
quark mass m c , and the NRQCD matrix elements. We take the QCD coupling constant to be a s = 0.21, which 
corresponds to the renormalization scale 5.3 GeV. For m Cl we use the next-to-leading order pole mass, which we take 
to be m c = 1.4 ± 0.2 GeV. The NRQCD matrix elements can be determined from the electromagnetic annihilation 
decay rates of J/ip and ip(2S). We use the values from the analysis of Ref. Q. For m c = 1.4 GeV, they are 
(Oi)is = 0.335 ± 0.024 GeV 3 and {0 1 ) 2 s = 0.139 ± 0.010 GeV 3 , where the error bars are those associated with 
the experimental uncertainties only. For other values of m c , the NRQCD matrix elements should be multiplied by 
(77i c /1.4GeV) 2 . 



TABLE I: Cross sections in fb for e + e~ annihilation into two S-wave charmonium states Hi + H2 at -Ebcam = 5.3 GeV. The 
errors are only those from variations in the pole mass m c = 1.4 ± 0.2 GeV. There are additional large errors associated with 
perturbative-QCD and relativistic corrections, as described in the text. 



H1 + H2 


a (fb) 


J/lp + J/lp 


6.65 ± 3.02 


J/i> + ij;(2S) 


5.52 ± 2.50 


^(2S)+V(2S) 


1.15 ± 0.52 


J/lp + Vc 


3.78 ± 1.26 


J/ip + Vc (2S) 


1.57 ± 0.52 


ip(2S) + Vc 


1.57 ± 0.52 


tP{2S)+ Vc (2S) 


0.65 ± 0.22 


rjc + Vc 


(1.83 ± 0.10)xl0 -3 


rjc + Vc(2S) 


(1.52 ± 0.08) xlO -3 


rj c (2S) + rj c (2S) 


(0.31 ± 0.02) xl0~ 3 



Our predictions for the double-charmonium cross sections for S'-wave states are given in Table [I] The error bars 
are those associated with the uncertainty in the pole mass m c only. The small error bars for rj c + rj c in Table [fl are 
a consequence of the value of m c being fortuitously close to a zero in the derivative of the cross section with respect 
to m c . The cross section for rj c + rj c is about 3 orders of magnitude smaller than that for J/ip + rj c . This suppression 



4 



comes primarily from the coupling-constant factor o? /a 2 . The cross section for J /ip + J/ip is larger than that for 
J/ip + r/ c . The suppression factor of a 2 /a 2 is more than compensated by the kinematic enhancement factor that 
scales as r~ 6 . The cross section for J/ip + J/ip is dominated by the photon- fragmentation diagrams in Figs. Ufa) 
and^b). For m c — 1-4 GeV, they contribute about 115% of the cross section. In Ref. 2], it was pointed out 
that there may be large perturbative-QCD and relativistic corrections to the production cross section for J/ip + r] c . 
There may also be large perturbative-QCD and relativistic corrections to the cross section for J/ip + J/ip. These 
corrections would affect not only the absolute cross sections in Table [I] but also the ratios of cross sections. For the 
J/ip + J/ip, J/ip + ip(2S), and ip(2S) + ip{2S) cross sections, the dominant photon-fragmentation contributions receive 
perturbative-QCD correction factors of about 0.39 and relativistic correction factors of about 0.78, 0.62, and 0.49, 
respectively |8j. 

The angular distributions da/d\x\ for m c = 1.4 GeV are shown in Fig. [21 At x = 0, the differential cross section for 
J/i/j + J/ip (normalized as in Fig.[5J) is smaller than that for J/ip + r] c by about a factor 0.66. However, the differential 
cross section for J/ip + J/tp is strongly peaked near the beam direction at x = 0.994, where it is larger than that for 
J/ip + r\ c by about a factor 9.3. The reason for the sharp peak is evident from the denominator in Eq. (J2J. 




FIG. 2: Differential cross sections da/d\x\ for e + e~ annihilation into J/tp + J/ip and J/ip + -q c at iSbcam = 5.3 GeV. The 
areas under the curves are the integrated cross sections 6.65 fb and 3.78 fb. There are large errors associated with QCD and 
relativistic corrections, as described in the text. 

The Belle Collaboration has recently measured the cross section for J/ip + r\ c by observing a peak in the momentum 
spectrum of the J/tp that corresponds to the 2-body final state J/ip + rj c |l| . The measured cross section is 

a[J/ip + r/ c ] x B[> 4] = (33±g ± 9) fb, (6) 

where B[> 4] is the branching fraction for the r] c to decay into at least 4 charged particles. Since B[> 4] < 1, the 
right side of Eq. © is a lower bound on the cross section to produce J/ip + rj c . This lower bound is about an order 
of magnitude larger than the predictions of NRQCD in the nonrelativistic limit 0,O|. Large relativistic corrections 
may account for part of the discrepancy ■ There may also be large nonperturbative corrections associated with the 
double-charmonium cross section being a significant fraction of the total cc + cc cross section . However the large 
discrepancy between the NRQCD predictions and the Belle measurement is still disturbing. This discrepancy can be 
decreased by taking into account the process e + e~ — » J/P , + J/ip- In the Belle fit. to the J/ip momentum distribution, 
the full width at half maximum of the r] c peak is about 0.11 GeV. Since the mass difference between the J/ip and the 
i] c is about 0.12 GeV, there are probably J/ip + J/ip events that contribute to the J/ip + r\ c signal observed by Belle. 

The Belle Collaboration also saw evidence for J/ip + Xco(l-P) and J/ip + r] c (2S) events. The i] c (2S) was recently 
discovered by the Belle Collaboration at a mass M^ c ( 2 s) = 3654 ± 6 ± 8 MeV |9j. Since the mass difference between 
the ip(2S) and the r/ c (2S) is only about 0.03 GeV, it is likely that any signal for J/ip + ri c (2S) in the J/ip momentum 
spectrum is contaminated by J/ip + ip(2S) events. A 3-peak fit to the data for the momentum spectrum of the 
J /ip gives approximately 67, 39, and 42 events with an accompanying r\ c ^ Xco(lP), or r] c (2S), with an uncertainty of 
12-15 events for each final state. The predictions in Ref. 0] for the relative cross sections for production of a J/tp 
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with an accompanying r/ c , XcoO-P)> or Vci^S) are 1.00, 0.63, 0.42, respectively. The observed proportion of events is 
compatible with the NRQCD predictions. 

There are also some unresolved puzzles in the inclusive J/ip cross sections that have been measured by the Belle 
and BABAR Collaborations [ToL Hl| . There are significant discrepancies between their measurements of both the 
production-angle distributions of the J/ip and the J/~0 polarization. There are also significant discrepancies between 
the measurements and predictions based on NRQCD [13 . These discrepancies may be decreased by taking into account 
photon-fragmentation contributions to the process e + e~ — * J/ip + qq, which allow the order-a 4 QED contribution 
to compete with the order-a 2 a 2 QCD contribution. The two experiments impose different cuts to decrease the 
background from the two-photon processes e + e~ — > e + e _ 7*7* and from initial-state-radiation processes such as 
e + e~ — > i/j(2S) + 7 and e + e~ — > J/ip + . These cuts remove substantial parts of the photon-fragmentation 
contribution to J/ip + qq. The differences between the experiments could be due to differences between the cuts. 
The discrepancies with the NRQCD predictions could be due partly to the neglect of the order-a 4 QED terms in the 
theoretical predictions. 

Another puzzling result comes from the measurement by the Belle Collaboration of the fraction of the inclusive J / ip 
cross section that is attributable to the final state J/ip + cc . That fraction is much larger than predictions based on 
NRQCD. The discrepancies between the predictions and the measurements would be decreased by taking into account 
the contributions to the processes e + e~ — > J/ip + qq and e + e~ — > J/ip + cc from photon fragmentation into J/ip. The 
experimental cuts allow the photon-fragmentation contribution to J/ip + cc to survive, but they remove substantial 
parts of the photon- fragmentation contribution to J/ip + qq. Thus, the net effect of the photon- fragmentation processes 
is to increase the ratio of J/ip + cc events to J/ip + X events. 

In summary, we have calculated the cross section for e + e _ annihilation into J/ip + J/ip. The calculated cross section 
is larger than that for e + e~ — > J/ip + r\ c by about a factor of 1.8 because it receives contributions from the process 
in which two virtual photons fragment into two cc pairs. The inclusion of this process in the analysis may decrease 
the large discrepancy between the Belle measurement of the production cross section for J/ip + r/ c and the predictions 
based on NRQCD. 
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In Ref. 0, we calculated the cross sections for exclusive production of two S-w&ve charmonium states in e + e~ 
annihilation. We made an error in the calculation by omitting a relative minus sign between diagrams that differ by 
the interchange of identical fermions in the final state. A similar error was made in the calculations in Ref. 0]. In 
\ what follows, we use results for the J/ip + r\ c , J/ip + -q c {2S), and J/ip + Xco(lP) cross sections in which this error in 
■ Ref. has been corrected. We show corrections in boldface type, except in display equations. 

The prediction for the ratio of the cross sections for J/ip + J/tp and J/ip + r/ c is given in both the abstract and the 
summary paragraph of Ref. jj. The correction of the error changes the prediction for the ratio from about 3.7 to 
about 1.8. 

In the expression for the differential ratio dR/dx for J/tp + J/ip, Eq. (3) should be corrected as follows: 

F(x) = r 4 a; 2 (l-a; 2 )[2 + 3r 4 -r 6 +4x 2 (l-r 4 )] 2 + 2(l + a; 2 )(l-a; 2 ) 
x [6 - 7r 2 + 4r 4 - r 6 + Ax 2 r 2 {\ - r 2 )} 2 + 2r 4 x 2 (l - x 2 ) 



x [3 - 4r 2 + 4r 4 - r 6 + 4a; 2 r 2 (l - r 2 )] 2 + 2r 2 [(2 - r 2 ) 2 (3 - 2r 2 + r 4 ) 2 
-x 2 (72 - 228r 2 + 359r 4 - 292r 6 + 130r 8 - 32r 10 + 3r 12 ) 
+4x 4 (9 - 60r 2 + 122r 4 - 112r 6 + 58r 8 - 14r 10 + r 12 ) 

+16x 6 r 2 (l - r 2 )(6 - llr 2 + llr 4 - 2r 6 ) + 64x 8 r 4 (l - r 2 ) 2 ]. (3) 



(N 

(N 
O 

Eq. (4) should be deleted. 

The expression for the differential ratio dR/dx for J/ip + rj c in Eq. (6) should be corrected as follows: 

<d : 

^m + Vo] = ^^(0 1 ) 2 5 [3a s r 2 + a (3 + r 2 )]V(l-r 2 ) 3 / 2 (l +a; 2 ). (6) 

•rH ! c 

| The QED terms increase the cross section by about 29% for i?bcam = 5.3 GeV. 

The error in the calculations produced errors in the numerical predictions for some of the cross sections given in 
Table I of Ref. 0. That table is reproduced here, with the corrected values shown in boldface type. The corrections 
decrease the " 3 Si + 3 Si" cross sections by about 24% and increase the " 3 Si + 1 5'o" cross sections by about 63%. 

The error in the calculations produced errors in Fig. 2 of Ref. pj. The corrected Fig. 2 is shown below. At x = 0, 
the differential cross section for J/ip + J/ip (normalized as in Fig. 2) is smaller than that for J/ip + rj c by about a 
factor 0.66. The differential cross section for J/ip + J/ip is strongly peaked near the beam direction at x = 0.994, 
where it is larger than that for J/ip + r\ c by about a factor 9.3. 

There are several additional corrections that must be made in the text. 

• The cross section for J/ip + J/ip is dominated by the photon-fragmentation diagrams in Figs. 1(a) and 1(b). For 
m c = 1.4 GeV, they contribute about 115% of the cross section. 

• The predictions for the relative cross sections for production of a J/ip with an accompanying r/ c , Xco{^P), or 
r) c (2S) are 1.00, 0.63, 0.42, respectively. The observed proportion of events is compatible with the NRQCD 
predictions. 

• The following sentence should be deleted: "If significant fractions of the J/ip + r\ c and J/ip + i] c (2S) signals 
are actually J/ip + J/ip and J/ijj + ip{2S) events, then the data would be more compatible with the NRQCD 
predictions." 
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FIG. 2: Differential cross sections da/d\x\ for e + e~ annihilation into J/ip + J/ip and J/ip + rj c at £bcam = 5.3 GeV. The 
areas under the curves are the integrated cross sections 6.65 fb and 3.78 fb. There are large errors associated with QCD and 
relativistic corrections, as described in the text. 
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TABLE I: Cross sections in fb for e + e~ annihilation into two 5*-wave charmonium states Hi + H2 at -Ebcam = 5.3 GeV. The 
errors are only those from variations in the pole mass m c = 1.4 ± 0.2 GeV. There are additional large errors associated with 
perturbative-QCD and relativistic corrections, as described in the text. 



Hl + H 2 


a (fb) 


J/V> + J/ip 


6.65 ± 3.02 


J/ip + iP{2S) 


5.52 ± 2.50 


tp(2S)+tp(2S) 


1.15 ± 0.52 


J/ip + rj c 


3.78 ± 1.26 


J/iP + Vc {2S) 


1.57 ± 0.52 


ip{2S) + Vc 


1.57 ± 0.52 


ip(2S) + Vc(2S) 


0.65 ± 0.22 


rjc + rjc 


(1.83 ± 0.10)xl0 -3 


rjc + rjc(2S) 


(1.52 ± 0.08) xl0~ 3 


rj c (2S) + ij c (2S) 


(0.31 ± 0.02) xlO -3 



